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Abstract A multiscale time integrator Fourier pseudospectral (MTI-FP) method is 
proposed and analyzed for solving the Klein-Gordon-Schrodinger (KGS) equations 
in the nonrelativistic limit regime with a dimensionless parameter 0 < e < 1 which is 
inversely proportional to the speed of light. In fact, the solution to the KGS equations 
propagates waves with wavelength at 0{e^) and (9(1) in time and space, respectively, 
when 0 < e <C 1, which brings significantly numerical burdens in practical compu¬ 
tation. The MTI-FP method is designed by adapting a multiscale decomposition by 
frequency to the solution at each time step and applying the Fourier pseudospectral 
discretization and exponential wave integrators for spatial and temporal derivatives, 
respectively. We rigorously establish two independent error bounds for the MTI-FP at 
0{t^/ h'^) and 0{e^ + h'"°) for e € (0,1] with t time step size, h mesh size and 
mo > 4 an integer depending on the regularity of the solution, which imply that the 
MTI-FP converges uniformly and optimally in space with exponential convergence 
rate if the solution is smooth, and uniformly in time with linear convergence rate at 
0(t) for e £ (0,1] and optimally with quadratic convergence rate at (9 (t^) in the 
regime when either e = (9(1) or 0 < e < T. Thus the meshing strategy requirement 
(or e-scalability) of the MTI-FP is T = (9(1) and h — (9(1) for 0 < e <C 1, which 
is significantly better than classical methods. Numerical results demonstrate that our 
error bounds are optimal and sharp. Finally, the MTI-FP method is applied to study 
numerically convergence rates of the KGS equations to the limiting models in the 
nonrelativistic limit regime. 
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1 Introduction 

Consider the Klein-Gordon-Schrodinger (KGS) equations in li-dimensions (d = 3,2,1) 
[2,22,42]: 


ihdt\f/{x,t) + - — A\f/{xd) + g^{xd)\j/{x,t) = 0, X e II 
2m I 

1 

^ar,0(x,f)-40(x,f)-f-^0(x,f)-g|v/(x,f)|2 = O, 


( 1 . 1 a) 

(1.1b) 


which represent a classical model for describing the dynamics of a complex-valued 
scalar nucleon field \j/ := i/(x,f) interacting with a neutral real-valued scalar meson 
field (j) := 0 (x, f) through the Yukawa coupling with 0 7 ^ g € K. the coupling constant. 
Here t is time, x £ is the spatial coordinate, h is the Planck constant, c is speed of 
light, mi > 0 is the mass of a nucleon and m 2 > 0 is the mass of a meson. 

In order to scale the KGS (1.1), we introduce 


7 , x=-, \fr{x,t)=x^J^\i/{x,t), ^(x,f) = 

ts Xs 


0 .? 


( 1 . 2 ) 


where Xs, ts and 0 * are the dimensionless length unit, time unit and meson field unit, 
respectively, satisfying ts = —^ and 0 ^ = with v = ^ = being the wave 
speed. Plugging (1.2) into (1.1), after a simple computation and then removing all', 
we obtain the following dimensionless KGS equations in c/-dimensions {d = 3,2,1): 


/5,i/(x,f) -f zii/(x,f) -I- A0(x,f)v/(x,f) = 0, X e : 




e-^40(x,f)-40(x,f)-t-^0(x,f)-A|v/(x,f)| = 0 , 


(1.3a) 

(1.3b) 


where e is a dimensionless parameter inversely proportional to the speed of light 
given by 

V Pi 

< 1, (1.4) 


„ V Xs h 

0<e:=- = — = -- 

c tsC 2cm\Xs 


are two dimensionless constants which are 


and ^ = 1 ;^ > 0 and A = e I 

independent of £. 

We remark here that if one chooses the dimensionless length unit Xj = 7 ^, ts = 

t = 27^ e = 1 in (1.4) and Eqs. (1.3) with 

e = 1 take the form often appearing in the literature [2,22,42]. This choice of x* is 
appropriate when the wave speed is at the same order of the speed of light. However, 
when the wave speed is much smaller than the speed of light, a different choice of 
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Xs is more appropriate. Note that the choice of x* determines the observation scale of 
the time evolution of the particles and decides: (i) which phenomena are ‘visible’ by 
asymptotic analysis, and (ii) which phenomena can be resolved by discretization by 
specified spatial/temporal grids. In fact, there are two important parameter regimes: 
One is when £ = 1 x, = ts = 5^ and 0, = then Eqs. (1.3) 

describe the case that wave speed is at the same order of the speed of light; the other 
one is when 0 < £ ^ 1, then Eqs. (1.3) are in the nonrelativistic limit regime. 

To study the dynamics of the KGS (1.3), the initial data is usually given as 


V/(x,0) = v/b(x), 0(x,O) = 0o(x), 5r0(x,O) = ^0i(x) 


X € . 


(1.5) 


where the complex-valued function y/o and the real-valued functions 0o and are 
independent of £. The KGS equations (1.3) are dispersive and time symmetric. They 
conserve the mass of the nucleon field 


\\¥{-,t)\\l 2 -= / lr(x,f)P^lx= f \\i/{x,0)\^dx= / |v/o(x)| 
Jm'' Jb'' Jmi 

and the Hamiltonian or total energy 

£-|o>,<^.|2 + |V0|2 + ^ 


E{t) := f 
./B‘ 

Iw' 


1 

£2 

= £(0), f>0. 


l<?>iP + |V0o|' 




i 0 r 


I'^'or 


|Vv/| 2 -A|v /|20 

-|Vv/op-A|v/op 0 o 


dx, f > 0, (1.6) 


dx 


dx 


(1.7) 


Eor the KGS equations (1.3) with £ = 1, i.e. 0(l)-speed of light regime, there 
are extensive analytical and numerical results in the literatures. Eor the existence 
and uniqueness as well as regularity, we refer to [3, 22-27,42, 45] and references 
therein. Eor the numerical methods and comparison such as the finite difference time 
domain (EDTD) methods and Crank-Nicolson Eourier pseudospectral method, we re¬ 
fer to [13,46,50,51] and references therein. However, for the KGS equations (1.3) 
with 0 < £ ^ 1, i.e. nonrelativistic limit regime (or the scaled speed of light goes to 
infinity), the analysis and efficient computation of the KGS equations (1.3) are math¬ 
ematically and numerically rather complicated issues. The main difficulty is due to 
that the solution is highly oscillatory in time and the corresponding energy functional 
E{t) = 0{e^^) in (1.7) becomes unbounded when £ —> 0. 

Eormally, in the nonrelativistic limit regime, i.e. 0 < £ <C 1, similarly to the anal¬ 
ysis of the nonrelativistic limit of the Klein-Gordon (KG) equation [38-41], taking 
the ansatz 

0(x,f)=e'>'/"VxT)+e-'^'/®'z(x,f)+o(£), xgR^ f > 0, (1.8) 

where z denotes the complex conjugate of a complex-valued function z, and plugging 
it into (1.3) and (1.5), we obtain the Schrodinger equations with wave operator as a 
limiting model 

id,yf{x,t) + Ayf{x,t) = 0, x € 

2ifidtz{x,t) + e^dttz{x,t) - Az{x,t) = 0, 


t > 0, 


(1.9a) 

(1.9b) 
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X 





Fig. 1.1: Plots of the solution of the KGS (1.3)-(1.5) with c/ = 1 for different e, where 
Im{t//} denotes the imaginary part of y/. 


with the well-prepared initial data [6,7] 


4x,o) = i 


0 o(x)--0 i(x) 

H- 


, dtz{x,Q) = -—Az{x,Q), v/(x,0) = i/o(x). (1.10) 
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In addition, after dropping the second term in ( 1.9b), formally we get the Schrodinger 
equations as another limiting model 

idt\ir{x,t) + Ayif{x,t) =0, (1.11a) 

2ijj.d,z{x,t) —Az(x,t) =0, xsK^, f > 0, (1.11b) 


with the initial data 


2(X,0) = - 


^o(x)--0i(x) 


V/(x,0) = i/o(x), xgW^. (1.12) 


This formally suggests that the solution of (1.3) propagates highly oscillatory waves 
with amplitude at (9(1) and wavelength at 0{e^) and (9(1) in time and space, respec¬ 
tively, when 0 < e <C 1. To illustrate this. Fig. 1.1 shows the solution of the KGS 
equations (1.3)-(1.5) with d = \, jx = X = \, \j/Q(x) = -^sech(x^), <l)o{x) = 
and 01 (x) = * for different 0 < e < 1 . 

The highly temporal oscillatory nature of the solution to the KGS equations (1.3) 
causes severe burdens in practical computation, making the numerical approximation 
extremely challenging and costly in the regime of 0 < e ^ 1. Recently, different nu¬ 
merical methods have been proposed and/or analyzed for the nonlinear Klein-Gordon 
equation in the nonrelativistic limit regime in which the solution shares similar oscil¬ 
latory behavior as that of the KGS equations (1.3), including the FDTD methods [9], 
exponential wave integrator Fourier pseudospectral (EWI-FP) method [9, 11 ], asymp¬ 
totic preserving (AP) method [21], stroboscopic average method (SAM) [18, 19] and 
multiscale time integrator Fourier pseudospectral (MTI-FP) method [ 8 ], etc. Among 
them, the SAM and MTI-FP methods are uniformly convergent (UA) for e € (0,1], 
while the FDTD, EWTFP and AP methods are not. Prom the practical computation 
point of view, the MTI-PP is much simpler and thus more efficient than the SAM 
method. The main aim of this paper is to propose and analyze a MTI-PP method 
for the KGS equations (1.3) in the nonrelativistic limit regime by adapting a multi¬ 
scale decomposition by frequency to the solution at each time step and applying the 
Pourier pseudospectral discretization and exponential wave integrators for spatial and 
temporal derivatives, respectively. Two independent error bounds will be established 
for the MTI-PP, which imply that the MTI-PP converges uniformly and/or optimally 
for 0 < e < 1. The MTI-PP method is also applied to study numerically convergence 
rates of the KGS equations (1.3) to its limiting models (1.9)-(1.10) and (1.1 1)-( 1.12). 

The paper is organized as follows. In section 2, we introduce a multiscale de¬ 
composition for the KGS equations (1.3) based on frequency. A MTI-PP method is 
proposed in section 3, and its rigorous error bounds in energy space are established in 
section 4. Numerical results are reported in section 5. Pinally, some conclusions are 
drawn in section 6 . Throughout this paper, we adopt the standard Sobolev spaces [1] 
and use the notation A < B to represent that there exists a generic constant C > 0, 
which is independent of time step T (or «), mesh size h and e, such that |A| < CB. 
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2 Multiscale decomposition 

Let T = Zif > 0 be the time step size, and denote time steps by = nr for n = Q,l,... . 
In this section, we present a multiscale decomposition for the solution of (1.3) on the 
time interval with given initial data at f = f„ as 


1 


0 (x,f„) = 00 W =0(1), a,0(x,f„) = ^0i"(x) , (2.1a) 

V/(x,f„) = vzo(x) = 0(1). (2.1b) 

Similarly to the analytical study of the nonlinear Klein-Gordon equation in the non- 
relativistic limit regime in [8,40], we take an ansatz to the variable 0(x,f) : = 0(x,f„ + 
s) of (1.3b) on the time interval [f„,fn+i] with (2.1) as [ 8 ,12] 


0(x,f„+i) =e'>*/®V(x,^)+e-'>*/®'i^(x,^)+r"(x,^), xs: 
Differentiating (2.2) with respect to s, we have 


0<s<x. (2.2) 




dsz’'{x,s)+‘-^z"{x,s) 


+ d,f^{x,s) 




Ifi. 


dsZ"{x,s) - -^Z’'(x,s) 


, xeR'', 0<i<T. (2.3) 


Plugging (2.2) into (1.3), we get forx G R'^, 0 < s < T and t//(x,f„ + s) =: i^"(x,s) 

e'W^' [e^d,,z"ix,s)+2i^d,z"{x,s)-Az"{x,s)] 

_l_e-'We^ [e^5siZ"(x,s) _2/^5iZ"(x,s) _z\z''(x,s)] 

+£2a,,,r"(x,.) - A^{x,s) + ^r"(x,^) = A | V'"(x,^)|" • 


Substituting (2.2) into (1.3a) and multiplying the above equation by e and 

gi/ji/e ^ respectively, we can re-formulate the KGS equations (1.3) for i/r" ;= \j/'{x,s), 
z" := z'’{x,s) and r" := r"(x,i) as 


/a,V/" + 4V'" + A (^e^z>"+e^^i^v^" + r>"^ =0, 

< lijidsz" + e^dssz" - Af =0, xgR'', 0<i<T, (2.4) 

e2a,,r"-4r" + ^r"=A|v/"p. 

In order to find proper initial conditions for the system (2.4), setting s = 0 in (2.2) 
and (2.3), noticing (2.1a), we obtain 

r z"(x,O)+i^(x,O) + r"(x,O) = 0o"(x), xgR'', (2.5) 

\ ^ [z'*(x,0)-?^(x,0)] +a,z"(x,0) + 5,F(x,0)+a,r”(x,0) = ®. 
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Now we decompose the above initial data so as to: (i) equate O (9(1) terms 

in the second equation of (2.5), respectively, and (ii) be well-prepared for the sec¬ 
ond equation in (2.4) when 0 < e ^ 1, i.e. 5iZ"(x,0) is determined from the second 
equation in (2.4), by setting e = 0 and i = 0 [6-8]: 


r z”(x,O)-f?'(x,O) = 0^'(x), [z"(x,0)-?'(x,0)] =0"(x), 

< 2/Ma,z"(x,0)-4z"(x,0)=0, xgR^, (2.6) 

[r"(x,0) =0, 5sr"(x,0)-|-(9,.z”(x,0)4-<3si"(x,0) =0. 


Solving (2.6) and noticing (2. lb), we get the initial data for (2.4) as 




(.^'(x)--<].f(x) 

M 


, a,z"(x,0) = -—4z"(x,0), xgK'', 


V/''(x,0) = i/o (x)i '■"(x,0) = 0, dy’{x,0) = -d,z'\x,0) - d,z’'ix,0). 


The above decomposition can be called as multiscale decomposition by frequency 
(MDF). In fact, it can also be regarded as to decompose slow waves at £^-wavelength 
and fast waves at other wavelengths, thus it can also be called as fast-slow frequency 
decomposition. After solving the decomposed system (2.4) with the initial data (2.7), 
we get y/"{x, t), z"(x, t), ds^{x, t), r"(x, t) and dsr''{x, t). Then we can reconstruct 
the solution to the KGS equations (1.3) at f = by setting s = T in (2.2) and (2.3), 
i.e. 

( r(x,Wi) = r"(x,T) =: i/o^^lx), 

I <^)(x,f„+i) =e'>^/^V(x,T)+e-'>V*^V(x,T) + r«(x,T) =: 0"+'(x), ( 2 . 8 ) 

|^5,0(x,f„+i) = ■^^"+'(x), xgR"', 

with 

^J'+'(x) [e^5,.z"(x,T)-|-/^z"(x,T)] [e^5,.i"(x,T)-//ii"(x,T)] 

+ e^dsr"{x,z). 


In summary, the MDF proceeds as a decomposition-solution-reconstruction flow at 
each time interval, and this makes it essentially different from the classical modulated 
Fourier expansion [14-17, 34-36] which only carries out the decomposition at the 
initial time f = 0. 


3 A MTI-FP method 

For the simplicity of notations and without loss of generality, we take /r = A = 1 in 
(1.3) and present our numerical method in one space dimension (ID). Generalizations 
to higher dimensions are straightforward and the results remain valid. We truncate the 
whole-space problem (1.3) into a finite interval Q = {a,b) with periodic boundary 
conditions, where due to the fast decay of the solution at far field, the truncation error 
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can be negligible by choosing a^b sufficiently large. In ID, the problem (1.3) with 
jj, = X = 1 collapses to 

idi\ir{x,t) + dxx'^f{x,t) + (l){x,t)\lf{x,t) =0, xG£2, f > 0, 

e^d„<j){x,t) - dxx(j}(x,t) + ^0(x,f) = \\l/{x,t)\^, xGQ, t >0, 

< ^ 

(j){a,t) = (j){b,t), d^(t>(a,t) = dj,(t>(b,t), \^{a,t) = \^{b,t) f > 0, (3.1) 

(j){x,0) = (j)o{x), d,(j){x,0) = v^(x,0) = v^o(x), x€Q = [a,b]. 

We remark that the boundary conditions considered here are inspired by the inherent 
physical nature of the system and they have been widely used in the literatures for the 
simulation of the KGS equations [10, 13,50,51]. 

Consequently, for ii > 0, the decomposed system MDF (2.4) in ID collapses to 

' id.xj/'’ + + =0. 

< 2ids^^ + e^dss^-d„f =0, a<x<b,Q<s<T, (3.2) 

^ e^d,,r"-d„r" + ^/‘=\r\\ 

The initial and boundary conditions for the above system are 

' z^ia^s) =z"{b,s), dj,z''{a,s) =dj,f{b,s), 

r^{a,s)=r^[b,s), d./'ia.s) = d,r^{b,s), 

\l/"{a,s) = \l/'\b,s), dx\i/''{a,s) = dxY\b,s), 0 < s < z; 

' z''(x,0) 1 [0o"(x) - (x)], 5,z"(x,0) = -la^.z''(x,0), 

r"(x,0) = 0, (3jr"(x,0) = -d^z"{x,0) - dsZ'‘{x,0), 

V^"(x,0) = i/o(x), a<x<b. 


In what follows, we present a numerical integrator Fourier pseudospectral discretiza¬ 
tion for the MDF (3.2) with (3.3), which applies the Fourier pseudospectral discretiza¬ 
tion to spatial derivatives followed by using some proper exponential wave integrators 
(EWI) for temporal discretizations in phase (Fourier) space. 

Choose the mesh size h := Ax = (b — a)/N with N a positive integer and denote 
grid points as Xj : = a + jh for j = 0,1,... ,N. Define 


XAr:=span<'e'^'(^ : x€Q, ...,- 1 

b — a 2 2 


N-\ 


Yn:= {v= (vo,vi,...,Viv) G C^+' : vq = Viv} , ||v||p =hj^ \vj\ 

1=0 


For a periodic function v(x) on Q and a vector v G l)v, let Pn-L^ {£2) —>■ be the 
standard-projectionoperator, and Ifj :C{£2) ^Xfj or Yfj — Xn be the trigonometric 
interpolation operator [31,48], i.e. 

N/2-l N/2-1 

iPNv){x)= ^ v,e'^'(^-“), (/^v)(x)= Y. 

l=-N/2 l=-N/2 


(3.4) 
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where v/ and v/ are the Fourier and discrete Fourier transform coefficients of the 
periodic function v(jc) and vector v, respectively, defined as 


V, = 


b — aja 


N-\ 


L 

j=0 


-ip.,(xj-a) 


(3.5) 


Then a Fourier spectral method for discretizing (3.2) reads: 

Findzj^, := rj(,(x,s), := %{x,s) € Xa? for 0 < s < T, i.e. 

N/2-1 N/2-1 

L r"^{x,s)= (';^),(.)e'>'(--«), 

l=-N/2 1=-NI2 

N/2-1 

W'^ix,s)= Y 

l=-N/2 

such that 

r id,\j/'!:f + d^^\j/^ + e‘'/‘^^PN{4f\j/N)+^^‘'^^^PN{^¥N)+PN{rN¥N)=^^ 

J 2/5sz5(, + e^dssZN - dxxz’^ = 0, 0 <i<T, a<x<b, 

— dxxr'lj + | 


(3.6) 


(3.7) 


Substituting (3.6) into (3.7) and noticing the orthogonality of the Fourier basis, we 
get 


2;(4)/(^) + e^(^7v)7 W + Mf(4)7(^) =0, 


0 < i < T, 
N 


(3.8) 




e ('‘at); (^)+ M/ +72 ('■5v)/(^) = (IV'wP)/(^)> -y < ^ < y - 1- 


For each —N/2 < I < N/2 — 1, we can rewrite (3.8) by using the variation-of-constant 
formula as 


Izn)i (s) =aiis) iz"fj)i{0) + e^bi (s) (0), 


0 < S < T, 




sin(a)/5) 77 




«>'-<»>+X 


sin(a)/(s — 0)) 


£ 2 ®/ 


(|V/J^|2),(0)7f0, 


(3.9a) 

(3.9b) 


0 ^,( 5 ) + (3.9c) 

+ /e-+^' jJ* e'(^'^ ? ) ® {d)de + i£ (;^); (0) uf 0, 


where 


, , A+e“V - Are“V 

ai{s):= -—r—-, bi{s):=i 


A,+ - Xf 


e"(V-V)’ 


0 < S < T, 


aX = -^ i± 


2g2 


(Ol = ^ 


2 £ 2 . 


(3.10) 
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Differentiating (3.9a) and (3.9b) with respect to s, we obtain 

{z'}j),{s)=a,{s){z"^)i{0) + e^b'i{s){z"^)i{0), 0 < i < T, (3.11a) 

(;^);(.) = cos((n,.)@(0) + £ ~ {J^),ie)de, (3.11b) 

where 


g«A, _gta+ 


’'<*> = V-V ■ = e^(V-A,-) 


Taking s = T in (3.9) and (3.11), we immediately get 


0<i< T. 


(4)/('^) = a/(''^)(^5v)/(0) + e^^/(T)(z;(,),(0), 
(' 4 );(T) = fl;(T)(' 4 ),( 0 ) + eV,(T)('^);( 0 ). 


In order to approximate the definite integrals in (3.9b), (3.9c) and (3. 1 lb) with s = T, 
we adapt the Gautschi’s type quadrature [8, 12,28-30,33, 33] 

[/(o) + 0/(o)] de = ^^1^/(0) + (^ ~ ~ ^ /(0) 

except the last term in (3.9c) which is approximated via a combination of the Gautschi’s 
and Deuflhard’s quadrature as [12,20,36] 


^ [( 10 ,( 0 ) + 0 (M( 0 ) 


de 


(Ai),(t) + t(A2),(t) 


(3.13) 


where Ai(jc, 9) := rj(,(x, 9)\l/'j!/{x,0) and A 2 (x, 0) := rj(,(x, 9)ds\l/^{x,0) for 0 < 0 < T 
and thus Ai(x,0) = A2(x,0) = 0 since r](,(x,0) = 0. Thus (3.9b), (3.9c) and (3.11b) 
with s = Z can be approximated as 




('^),W ^ +P/(T)(j^) 00 )+^,(T)(j^^ 

(OOOt) ^«cos(tu,T)(7^);(0)+p;(T)(I^),(0)+^;(T)(I^);(0), (3.14) 

(rw)/(T) ~ e-'>^(v/)^),(0 )+c+(t)(z)(,v/;^),(0)+^/+(t)(z« v'w)/( 0 ) 


IT 


+G iT)iz'^%),{0)+d, {T)iz"^%)i{0) + -^iAi),{T) + T{A2),iT) 


where pi{t), ^/(t), p)(t), ^)(t), cf{T) and df{T) are given in Appendix A. For the 
derivatives, we can compute them as 


(]^)',( 0 ) = 2 (Re{v^v/;^}),( 0 ), + (z^;)),(o), 
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where Re{z} denotes the real part of a complex number z and for — f <Z<f-l, 

(^Mo) = Inm= +/(m)/(o), 


which are approximations of dsz''{x,Q) = —^dxxZ"{x,Q) and dsy/"{x,0) = 

+ /0"(a:,O)v/'"(x,O), respectively [8]. Inserting (3.12) and (3.14) into (3.6) with setting 
s = T, and noticing (2.8), we immediately obtain a multiscale time integrator Fourier 
spectral method based on the MDF (3.2) for the problem (3.1). 

In practice, the integrals for computing the Fourier transform coefficients in (3.5), 
(3.9) and (3.11) are usually approximated by the numerical quadratures [9, 11,48] 
given in (3.11). Let ^>", 4>'' and be approximations of (j)(xj,t„), dt(j)(xjjn) and 
\f/{xj,t„), respectively; and be approximations of z''{xj,z), 

dsz"{xj,z), r'^ixj, t) and t), respectively, for j =Q,1,...,N. Choosing (pj = 

0 o(jc;), 4’“ = <l)i{xj)/£^ and = y/o{xj) forO < j <N,a multiscale time integrator 
Fourier pseudospectral (MTI-FP) discretization for the problem (3.1) reads for n>Q, 


' ^."+1 = +e-'^/'="z"+i +R';+^ 7 = 0,1,... ,iV, 


^n+i ^ 2''+! + -LZ"+' 1 +e-'^/'' Z"+^ - -LZ' 


N 

T“ 


1 '''+' = ^ (•?«+!),e' 


c2 J 


>7+n Jfii{xj-a) 


N 

"T 


7 «+l 


p 2 J 


+ R 


in+l 


(3.15) 


where 


ylt+l 


yll+l 


y -1 
/- ^ 


l—Y 


R”+^^ ^ 

J- ^ 

/—Y 

N_l 

R"+^= £ 

J- ^ 

/—Y 


and for -N/2 < I < N/2 - 1, 


(Z«+i),=fl,(T)(Z0), + e2l,,(T)(Z0)„(Z«+i),= fl;(T)(Z0), + e2l,;(T)(Z0)„ (3.16a) 

(Oi 

{F^),=C0^[C0iT)(¥),+p\{T)(^),+2q\{T){Rt{^'¥-}),, 


(3.16b) 

(3.16c) 


CF-'+i), =e-'>/ ''('F"), +c,+ (T)(Z'>'f'«), +£/+(t) 


(Z 0 'F«), + (ZO'F"), 


+ cJ{x){Z^W-),+di{x) 


+ (Z“'f"')/ 






(3.16d) 
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with 



(3.17b) 


(3.17a) 


This MTI-FP method for the KGS equations (3.1) (or (1.3)) is explicit, accurate, 
easy to implement and very efficient due to the discrete fast Fourier transform. The 
memory cost is 0{N) and the computational cost per time step is 0{N\ogN). 

4 Uniform convergence of MTI-FP 

In this section, we establish an error bound for the MTI-FP method (3.15) of the 
KGS (3.1), which is uniformly for £ G (0,1]. Let 0 < T < T* with T* the maxi¬ 
mum existence time of the solution of the problem (3.1), motivated by the formal 
asymptotic results (1.8) and (1.9) as well as the numerical results (cf. Fig. 1.1), we 
make the following assumptions on the solution of the problem (3.1) - there exists 


an integer iuq > 4 such that <f) G {Q)), Xj/ G C{[0,T]-,H'p°^^(Q)) fl 

C* ([0, r] ;7 /™“ (72)) n ([0, r] ;(72)) and 



(4.1) 


where77;;’(72) = |/(x) g 77'"(72) | /«(a) =/W(l7),k = 0,1,... ,m - 11 C 77'"(72). 
From the first equation in (3.1), i.e. the Schrodinger equation, it is easy to see that 

idtp{x,t) + dx \\i7dx\l/{x,t) — \lfdx\j7(x,t)] =0, x G 72, f > 0, 
where p{x,t) := |i//(x,f)p. Then under the assumption (4.1), we have 


p €C{[0,T];H’"«+^{Q)) nC^ ([0,7’];777(72)) nC^ , 

11*^? P 1) 7: = 0,1,2. 


(4.2) 



DenoteC,|) —maxo<£<i |ll0lli”“([o,r];//'"o+''(r2)) j — 

niaXo<£<l III '/lli'”([o.7’];//"’0+2(i2)) }■ 



0",..., 0'j^) G Y]\i (n > 0) are the numerical solution obtained from the MTI-FP 
method (3.15)-(3.17), denote their interpolations as 


0;(x) := (7^0")(x), Wii^) ■■= := (7^0")(x), x G 72, (4.3) 


and define the error functions as 


e 4 ,{x) ■■= - 0/'(x), e"^{x) := \l/{x,t„) - Xj/'/ix), 


— T (4.4) 

e^x) := d,<j){x,tn) — ^j{x), xG72, 0<n<—; 


then we have the following error estimates for the MTI-FP method (3.15)-(3.17). 
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Theorem 4.1 (Error bounds of MTI-FP) Under the assumption (4.1), there exist 
two constants 0 < ho < I and 0 < To < 1 sufficiently small and independent of e such 
that for any 0 < £ < 1, when 0 < h < ho and 0 < T < To, we have 


^2 + \K\\h^~+^^ 


</i '"0 + £ 2 , 

//2 


|| 0/"||//2 < Q + 1 , II V' 7 II //2 < 1 ; 


T 

0 <n < —, 

llnll <G±1 


(4.5) 

(4.6) 

(4.7) 


Thus, by taking the minimum of the two error bounds (4.5) and (4.6) for £ G (0,1], 
we obtain a uniform error bound with respect to e € ( 0 , l]for 0 <n < j 


^2 + F'/ll //2 + e 


7/2 


< /i'"« + min I < ffi'o + t. (4.8) 
0 <e<l ^£2’ j 


In order to prove the above theorem, we introduce 

:= (iV0)(x,f„)-0;(x), e'^_^(x) :=fV(5,0)(x,f„)-07(2c), 

— T 

e'wNi^) ■= {PNV^){x,tn)-W'/{x), xen, 0<n<—. 

T 

Using the triangle inequality and noticing the assumption (4.1), we have 

lk '^||//2 < || 0 (x,f„) — {PN^)(x,t„)\\[ji + ||e 5_^||//2 < + ||e^;v||/f 2 , 

lk^/ll //2 < ||r( 2 c,f«)-(^VV^)(x,f „)||;^2 + ||e '^_^||^2 </i™'> + ||e^_^||^ 2 , 

/l™ 0+2 

||e ^||//2 < ||0(x,f„) - {PN^){x,tn)\\fj 2 + ||e^,Ar ||//2 < —^ + ||e^,A,||H 2 . 


(4.9) 


(4.10) 


Thus we need only obtain estimates for 11 e'^ ^ 11 ^ 2 , 11 e^ ^ 11 ^^2 and 11 e^ ^ 11 ^^ 2 , which will 
be done by introducing the error energy functional 


^ \y,N^,N 


:= £ 




H2 


d,el 


'xCAN 


1 


7/2 £2 




H2 ^ lkv,A'll772 j ’ 

(4.11) 

and establishing the following several lemmas. 

For any v G Fiv, denote v_i = v'Ar_i and v'a^+i = vq, and define the difference 
operators 5+v € Uv and 5^\ G Vjv as 


5 ^ V j = 


Vj+\ ■ 


2 _ V;+l-2Vj+Vj_l 


^xVi = 


h^ 


; = 0,i,...,iv, 


and the norms as 
|2 ||.,||2 


V||f 1 = ||v||f 2 + II 5>||f2, ||v||f 2 = ||v||f 2 + II 5>||f2 + ||57v||f2, (4.12) 


l|v||p =/j |V;f, 

7=0 


\^Ml = hY\^xyj\\ ll5Ml?2 = /7i:i5.^;i , 

7=0 7=0 


77-1 

2.,ll2 _ V- is2^,.|2 


with 
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then we can have the following estimates [7, 8 ] 

< ||v||yj < ||/Arv||^i, ||/ivv ||^2 < ||v||y ,2 < ||/ivv||^ 2 , Vv £ Fat. (4.13) 

Lemma 4.1 (A new formulation of MTI-FP) For the numerical solution 0”+^ and 
0”^' (n > 0) obtained from the MTI-FP method (3.15)-(3.17), we have 

N/l-l N/l-l 

^’l+\x)= Y. 0;'+‘W= Y (4.14) 

1=-NI2 1=-NI2 


where 


=cos(®,T)(^, + +p,(T)(j^), +2^KT)(Re{^"})/, 

(Oi 

(0«+i), =cos(®,t)( 0«),- ®/sin(®,T)(^,+ /7',(T)(|'f""|2),+2^;(T)(Re{^'P"}),. 
Proof For / = —— 1, from (3.15), we get 

(0^), =e?(Z^),+e^? +{R^)i, (4.15) 

(0"+i),=e? (Z"+i),+ ^(Z'’+i), +e^? +(«"+')„ 


and from (3.17), we obtain 


1 


(Z0), = - {<P"),-ieH<P"), , (R0), = -(Z0),- zo . 


2 L 

Noticing (3.10) and (3.4), for v G Yn, we have 


(4.16) 


N 


N N ~ 

a/(T) =a_,(T), bi{T)=b_i{T), 1 = --,...,--!, (v), = 


V-/, |/|<y-l, 


V/, / = - 


A 


Plugging (4.16) into (3.16a)-(3.16c), and (3.16) into (4.15), we get 
( 0^)1 =Re|e?a,(T)| (^, + e2lm|e?a,(T)| (0"); + e^e?Z 7 ,(T)(Z 0 ), 


+ e2e 7^bi{z)\^Z^j^ + 
+ 2 ^,(T)(Re{^'f"'})„ 


sin(w;T) rr 


(Oi 


{R0), + p,iT){\^'^\2), 


I 4- 


=e2lm|eE2aJ('r)+ -^ee^ai{T:)j {<P'')i+e^e 
+ Re i ee^a'fx) + ^e^a/(T) i (^’")/ + 


b',{T)--^bi{T) 


(Z^)/ 


+ cos(cu,T)(RO),+yKT)(l'R'f)z + V;(T)(Re{'f">'f"'})„ 


where Im{z} denotes the imaginary part of a complex number z. Then (4.14) can be 
obtained by combining the above equalities and (3.10). □ 
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For the solution (j)" (x, s) := 0 (x, + s), V^" (x, s) := t„ + s) of (3.1), we have 


£ (p” {t„ + s) + Hi ^i{t„ + s )-\—2 ~ s>0 (4.17a) 


iw'l{fn ~ +'S) + (0V^)/(fn + ■*) — Oj / € Z. 


(4.17b) 


By applying the variation-of-constant formula to (4.17) and noticing the ansatz (2.2) 
with p"(x,s) := I (/(x,?,, +s)p, we get 

coi Jo e^coi 

(4.18a) 

^'(f„+i) =cos(®/T)^'(f„)-cu,sin(a)/T)0/(f„)+ [ ^^\ p"),{9)d0, 

Jo £ 

(4.18b) 

Vi,(f„+i) =e^'^^V>/(f„) + ie^'^’' f (4.18c) 

Jo 

Jo Jo 

Combining the above results, we define local truncation error functions as 


N/2-\ iV/2-1 

^;(^)= L ^;(x)= i: 

I=-NI2 l=-N/2 

N/2-1 _ _. 

^;w= E ^ei 2 , 

l=-N/2 


(4.19) 


where 


:=0/(Wi) - [cos(®,T)^,(f„) + + p,(T)(|V^),(0) 


0 )/ 


+ 2 ^,(T)(Re{v/"r})/( 0 ) 


(4.20a) 


= -(0lSm{C0iT)piit„)+COs{C0iT)pl{t„)+p',{T){\\t/"\^)l{0) 


+V/(T)(Re{rr})/(o) 


(4.20b) 


-flf+(T) [(fr);(o) + (?r),(o)l -diiT) \^r^),io) + ^r)iio) 


- Y(r+V")/(o) - ^(r+>")/- 


(4.20c) 

Here we introduce three auxiliary functions as 
<P"W = E 3"(x) = ^ (yO,e'>'(^-«), 7"+'(x) = ^ 
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with 

(^")/ =-^sin(^,2T)(v/'')/(O)+/(0"r")/(O), ( 3 ")/ =^sin(^/T)(z”),(0),(4.21) 

Similarly to the case of the KG equation [ 8 ] with details omitted here for brevity, 
we can prove a prior estimates for the solution of the problem (3.2) with (3.3). 

Lemma 4.2 (A prior estimate of MDF) Let z'’{x,s) and r’^{x,s) be the solution of 
the MDF (3.2) with initial boundary conditions (3.3). Under the assumption (4.1), 
there exists a constant Ti > 0 independent of £, h and 1 (orn), such that for 0 < T < Ti, 


11^ IIl“([0,t];//"' 0+2) + ll^s^ IIl“’([0.t];//'"0+‘) + IIl"“([0,t];//'"0) 






< e 


2-2k 


<1, (4.22) 

k = 0,1,2. (4.23) 

Based on this result, when 0 < T < Ti, it is easy to see that (p", 3 " G Hp'° {£2), y"^ * G 
Hp(^Q) and we have the prior estimates 

ll<P"ll/r4<l, ||3"ll/r4<l, 

With the above prior estimates, we can estimate the local truncation error functions. 

Lemma 4.3 (Estimates on local errors) Under the assumption (4.1), for any 0 < 
T < Ti with Ti given in Lemma 4.2, we have two independent bounds as 


< .j.2g2 


0< «<--!. (4.24) 

T 


IP (pn pn pn\ ^ ^ e (pn pn pn\ 

^ ^0 y ^ y 

Proof Subtracting (4.20) from (4.18) and by using the Taylor’s expansion, from 
(4.18a) and (4.20a), we get 

(1P,= j\l-s)(^)[{ds)ds\dd+2q,{x){B^\, (4.25) 

JO c CU/ ^«/0 


where 


Bo(x) = Re{ v/"(x,0)(5si/"(x,0) -(?")}, 
Noticing (A.l), we obtain 

r2 rx {-I 




(lrP)/( 0 * 


dsdd 


e^coi 


X £ £2. 


(Bo), 


Noting / < 1 for / = -f - 1 and(4.1), (4.2), we get 


E-COi 








< 


T^ + T^\\d,ri;0)-<P%2 


0 < T < Ti. 


(4.26) 
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By Parseval’s identity, we have 

Njl-X 2 

\\d,r{; 0 )-<P% 2 = L (i+^f+nt) i%r),{0)-(r), 


l=-N/2 


< L W - T - 

1=-NI2 ^ 


A'/2-l 


l(rMo)i"<T2 £ M/"i(r)/(o)i' 


\=-Nl2 


This together with (4.26), we get 

m\\l2<T\ 0 <T<T,. 

Similarly, noting ^ ^ for 1 ■ ■ ■, f — 1, we obtain 


(4.27) 


(4.28) 


\/l+rW " " 

t'’ . ^ T" 


||^''||^ 2 <—, o<T<Ti. 

For from (4.18c) and (4.20c) and by the Taylor’s expansion, we have 


_l_/e^'^/^ / e‘(^' [ (1 — s)(z"v/"), (0s)^f, 

Jo LJo 


(4.29) 


(4.30) 


dO 


+ dl{z) (Bi), + (B 2)/ +di {t:) {Bt,)i + {B4)i + —(B5), + y(Bg), 


/T^ '■ 

+ Tl 


+ (bT),] - ^(7^),(t)+ i^%'>?(®-^)(7^),(0)«f0, 


where 


Bi{x) = [dsz" (x, 0 ) - 3 " (x)] Xj/" (x, 0 ), B 2 (x) = z” (x, 0 ) [B,. \l/"{x,0) - (p" (x], 

B 3 (x) = [B,.i"(x,0)-^(x)] V^"(x,0), B 4 (x) =7'{x,Q)[dsy/'‘{x,Q)-(p'‘{x)], 

B 5 (x) =r"(x,T)Bi. 5 V^"(x, 0o), Be{x) = [r"(x,T) - 7 "+^(x)] V^"(x,0), 

Bj{x) = r"(x,T) [B,vr"(x,0) - (p"(x)], B 8 (x) = [r"(x,T) - 7 "+'(x)] (p"(x), 

with 00 € [0, t]. For the last two terms in (4.30), we have 

/J e'^''^(®^'')(7^),(0)£/0-y(7^),(T) 




B0. 
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Then by using the triangle inequality, we obtain 

rt /*1 




(z«r); ids) 


(z«r)/ (Os) 


dsdd 


+ r 


+ r 



+ r 


+ r 


L 


(Bi), + (B 4 )/ 

d^ 


+ T" 


(B5), +T (Be)/ 


^02 L' 




£/ 0 . 


Similarly as (4.27), we find from (4.21) 

\\Bsz" —3”||ef2 ^ (ll^llL“([o.r];//S) + ^^ll'^'^llL°°([o,r];//S)^ ^ T. 

Again, similarly as (4.27) and noticing (3.9b), we have 


lez 

By the assumption (4.1) and Lemma 4.2, we get 


sin(a)/T) 


(Oi 




+ + T^llf +1 - r"(., T)||^, + T® |k>"|l^.([ 0 .r];//«) 

+ T®||Bj('^’/")llL“”([O.T];//4) + T®||<9.5i(r'V")|lL=<.([o,T];//2) ^ (4-31) 

Plugging (4.28), (4.29) and (4.31) into (4.11) with and = 

we immediately get the first inequality in (4.24). 

Next we show how to get the second estimate in (4.24). By using Taylor’s expan¬ 
sion, truncating at the first order term and integrating by parts, we have 

^ r^ smicoiiT-9)) r ,1 _2^,(^)(Re{^«}),(0), 

^ ' Jo e^coi iJo J 

=‘}M [ (lr''P)/(T*)'^s-/ qi{ 0 )\[ {\¥''\^)i{ds)sds 
Jo 7o Ljo -I 

-^,(0)(|rf)'(0)-2^,(T)(Re{'r(j!)"})/(0). 

From (A.l) and noticing o'/(0) = 0 and k/('r)| < ze'^, we find 


d9 


U^\\h2 ^ T e [ll^'PllL-»([0,r],//2) + \\^‘rP\\L-{m.H^-) + II V^'P"llL-([0.r],//2)_ 

< 0<T<Ti, 0<e<l. (4.32) 

Similarly, we can get 

llkv^;il« 2 <TV and || 4 "||^ 2 <Tk 0<T<Ti. (4.33) 
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In order to estimate ^1, we introduce 


where 


with 


N/2-l N/2-1 

r(x)=52 J'’{x)=Y, 

l=-N/2 l=~N/2 


I" = Gh-dH^) [(3"r),(0) + (z”(p'')/( 0 ) 
Ji = Gh - di (T) [ (fr), (0) + (1V), (0) 


dO, 

de. 


(4.34) 


(4.35) 


Jo 

/■i —— 

/ {z''\tf")i{9s)ds 
LJO 

Jo 

/ {z''\tf")ii9s)ds 
iJo 


(4.36) 


Applying the Taylor’s expansion in (4.18c) and (4.20c) and truncating at the first 
order term, we obtain 

(r+V")/(0) + T(r+‘9'')7 • 
(4.37) 

From (4.36), integration by parts, we get 


Gl, = ie 


ie f {dsZ'’\l/"),{9s) + {z'’ds\l/")i{0s) 

Jo Jo L 

d^{z)J {dsZ"\l/")i{'Zs)ds — J d^{9)J {dsZ"\tf'')i{9s)sdsd9 

f f (z"d„nf")i{ 9 s) +{z”^’^\l/")i{ 9 s) 

Jo Jo ^ 


dsd9 




dsd9 


: d^ (t)^ + !(z”5„V^");(Ts) 

- d/{9)j^ {dsZ"\l/"),{9s) + i{z"d„\l/'’)i{9s) 


sdsd9-G'li, (4.38) 


where 


Gli=c 


[\f^^)A9s)dsd9. (4.39) 

Jo Jo 


Noticing the ansatz (2.2), we get 


G?./ = e 


+ 


Jo Lio 


J 9{\z"\^y/")i{9s)d9 + J 9(z"r"^f'')i{9s)d9 


ds. 
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Using integration by parts and the triangle inequality, we have 

Jo L JO 


((z«) 2 v/")/( 0 ^) 


d9 


+ £■ 


7 ’ 




d9+ [ 

Jo 


£/0 




(4.40) 


Combining (4.34)-(4.40) and noting |^^7('f’)l ^ Lemma 4.2, we get 


P1I//2 [||^s2>"IIl"“([0,t];//2) + lk"^«V'"llz.-([0.T];//2) + IIIIl“([0.t];//2) 

+ ||^iZ"^sV'"llL“([0,T];//2) + II^s^”^«V'"IIz."“([0,t];// 2) + lk”^sxc V^" |1 l">.([0,t];//2) 


+ Il(/7>"llz.==([0.r];/r2)_ 

+ T||z"r''v/"||i.([0_,].^2) 


+ Te-" 

< Te7 


I|3"'/"IIl"“([0.t];// 2) + lk"^"llL""([0.T];//2) 


(4.41) 


In order to estimate 7", denote Ne := when |1| < A^£, we have |/J; | < which 

indicates that (t)| < Te^. Similarly to the estimate of /", we have 


^ (i+Mf+M/)l777<TV. 

|;|<We 


(4.42) 


On the other hand, when |/J;| > we have (t)| < T^. Due to the regularity of 
the solution, i.e. the assumption (4.1) and Lemma 4.2, we can get 


^ (1+M?+Mf)|ir7 

l/U^E 


< L 

|/|>We 

< 


ff 

Jo Jo 


||57^>")lll==([o,.];/r4) + ||3>"(-,0)11^4 + ||z"(-,0)(p"||^4 


(z«r),( 0 ^) 


dsd 9 + 1 


(3«r)/(0) + (z'’<P”),(0) 


«Ii 2 


<tV 

< t2£7 


Il^.■(^"r”)lli^([0,.];/r4) + :3lk>"ll2^([0.r];//4) + lk>'>"ll2^([0.r];//4) 


(4.43) 


From (4.42) and (4.43), we obtain 

ll■^"ll //2 <Te^, 0<T<Ti, 0<e<L (4.44) 

Combining the above estimates and noting (4.37), we have 

+ T*||f+V"||J,:<lV. (4.45) 

Plugging (4.32), (4.33) and (4.45) into (4. 1 1), we get the second estimate. □ 
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Defining the errors from the nonlinear terms as 


N/2-\ N/l-l 

l=-N/2 l=-N/2 

N/2-1 

v;(x)--= E n>0, 

l=-N/2 


(4.46) 


where 


(j?P,=p/(t) [drf )/(o)-(|'f"f )/J+ 2 ^,(t) 
(^)rp'M [(iV^)z(o)-(mJ+2^;(T) 
(li)l =C;+(T) [(?r),( 0 ) - (Z^),] +C^(t) 


(Re{vr"(p"}),(0)- (Relf"*•?"}), 

(4.47a) 

(Re{^"}),(0)-(Re{^'f"'}), , 

(4.47b) 


(z'>«),(0)-(ZO'R«), 


(4.47c) 


+ (7,+ (t) 
+ c/, (t) 


(3"r)/(0) - (ZO-R")/ + (2"<P")/(0) - (ZO-R"), 

(fr)i(0) - (5“'R")/ + (1V)/(0) - (^'R"), 


!T 

+ 2 


(yn+ly/'fl)^( 0 ) _ ^ ~ (^"+'^"0/ ) 


then we have 


Lemma 4.4 (Estimates on t]^', 77 " and f/^ Under the assumption (4.1) and assume 
(4.7) holds (which will be proved by induction later), we have for any 0 < T < Ti, 

(^5v,(|u^5v,v/>^5v,i|i) "I-^ 2 —’ 0<n<—— 1. (4.48) 

Proof From (4.46) and (4.47a), we have 

< T||/7v|r"(-,0)|2-/^|'R"|2||^2+ Td|/7v(r(-,0)(p")-/7v(^'R")ll/72+ T/7'"« 

< Till V/'H^0)|2 - |'R«|2||y2 + T^lr (^0)(p"-■^'R"||y ,2 + T/r'"« 

<T||c;||y ,2 + Td|(p"-'R"||y .2 + T/7'”° 

<T||4||y,2 + T2||4V'''(^O)-4'R'd|r,2 + Td|0>"(-,O)-0"'R"||K.2 + T/r'"'> 

<T|k;ii,,2 + T||c;ii,,2 + Tdi4r(-,o)-ai,'R"ll,,2 + T/7-o 

<T|k ;||,,2 + T|k "^||,22 + T/7"-«, 


(4.49) 
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where the operator <9^ on a function /(x) is defined 


as 




Combining (4.49) and (4.10), we get 
Similarly, we can obtain 


(4.50) 


(4.51) 


Similarly, we can estimate all the terms of Tjj^ in (4.47c) except the last term. For the 
last term, we have 

^ ( 1 + 4 + 4 ) (r+V")/(o)-(/?”+''f"’)/ 

l=-N/2 

< ||/^(7”+V''(-,o))-/4/?"+V"))II4 

< \K\\l2 + \\f^^-lNR''+^\\l2+h^'’'°. (4.52) 

Noting (4.21) and (3.16b), we get 


< (z”),(o)-(z'>), +4 (t)|. {\r\^),{o)-m^)i 


+ k/WI' 




Combining the above two inequalities, we obtain 

||f+i-/^7;"+'l4<||4i4 + e2|l4l4 + ||4i4+/r2^^ 

Plugging (4.54) into (4.52), we get 

^ ( 1 + 4 + 4 ) (r+V"),(o)-(/?"+''f'")/ 

l=~N/2 

Similarly, we can obtain 


(4.53) 


(4.54) 


(4.55) 




N 12 
/=-A'/2+l 

<ll4llH2 + lk';ilH2 + e1^;il^2+/^""'“. 
Combining (4.55), (4.56) and (4.47c), we have 


(4.56) 


11 


;\W- < ^ {\\e;\\H2 + e^WnWn^ + Kh^+h’’’^) 


5, T 


{\\elM\\H2 + eHelM\\H2 + 114 ,^ 11 //^ • 


Plugging (4.50), (4.51) and (4.57) into (4.11), we get the estimate (4.48). 


(4.57) 

□ 
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Proof of Theorem 4.1. The proof will be proceeded by the energy method and the 
method of mathematical induction [6,7, 9,11, 12]. For « = 0, from the initial data in 
the MTI-FP method (3.15)-(3.17) and noticing the assumption (4.1), we have 

11411//^+11411//^+^'11411//^ 

= \\<Po-In<Po\\h^ + \\Vo-InWo\\h^ + 1101 -In<Pi\\h2 < < h'”°. 

In addition, using the triangle inequality, we know that there exists /q > 0 indepen¬ 
dent of e such that for 0 < h< hi 


110/ II //2 < || 0 oll //2 + ||e 0||//2 <C^ + l, 110/ II //2 < 
Il'Z/Ilz/a < II V4)||//2 + ||4II//2 < C^+ 1. 




Thus (4.5)-(4.7) are valid for n = 0. Now we assume that (4.5)-(4.7) are valid for 
0 <n <m—l < T jx—X. Adding (4.20) and (4.14), we have 


= cos(o,t)('4), + 


(4+*), = e-'>'"(4), + 

Using the Cauchy’s inequality, we obtain 


(4.58a) 

(4.58b) 

(4.58c) 




< (1 + t) 


cos(cu/t)(4)4 


+ ■ 


I + T 


(4 ), <(1 + t) cos(®/T)(4)^-®/sin(®/T)(4)^ 


2 l-t-T 




<(1+T) (4)/ 


2 l-f T 
H- 


(4)/+ ( 4 )/ 


(4)/+( 4 )/ 

(4.59a) 

(Ip/+ (4), ’ 

(4.59b) 

(4.59c) 


Multiplying (4.59a)-(4.59c) by e^{l+^f + ^l) and ^(1 

+ 11^), respectively, then adding them together and summing up for / = —N/2, 
... ,N/2 — 1, we obtain 


Using the Cauchy’s inequality, we get 

/ 2 >/ f?+l n+1 -n+lx /v/ n n •« \ 


1+T 




(4.60) 
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Inserting (4.48) and the second inequality in (4.24) into (4.60), we get 




^l^2mo 


(4.61) 


Summing the above inequality for 0 < n < m — I and then applying the discrete 
Gronwall’s inequality, we have 




< 


/j 2 «io 


6 ’ 


£ 

Similarly, by using the first inequality in (4.24), we obtain 

/, 2 mo 


n n ^ 


el 


(4.62) 


(4.63) 


Combining (4.11), (4.62) and (4.63), we get that (4.5) is valid for n = m, which 
implies 

\W;\\H2 + \\e';\\H2+e^\\e';\\H2<h"'<> + T. 

Using the triangle inequality, we obtain that these exist /12 > 0 and T 2 > 0 independent 
of e such that 


||0ni//2 < \\<j>{-,tm)\\H2 + < Q + 1, 

llr/'ll//2 < llr(->fm)||//2 + ||e“||222 <Cv, + l, 0<h<h2, 0<T<T2, 

ll0rL2 < ||5,0(-4.)L2 + |k;il2,2 < ^. 

Thus (4.7) is also valid for n = m. Then the proof is completed by choosing To = 
min{Ti,T 2 } and /iq = min{/ri,/r 2 }. □ 


Remark 4.1 Here we emphasize that Theorem 4.1 holds in two/three dimensions 
(2D/3D) and the above approach can be directly extended to 2D/3D without extra 
effort. The only thing needs to be taken care of is the Sobolev inequality used in 
Lemma 4.4 in 2D/3D 


ll“llL°°(r2) ^ Il“ll//2(i2), ll“llw*'P(i2) ^ Il“ll//2(i2)i 1 < p < 6, (4.64) 

where is a bounded domain in 2D/3D. By using the assumption (4.7), Lemma 4.4 
still holds in 2D/3D. Thus (4.7) and the error bounds can be proved by the method of 
mathematical induction since our scheme is explicit. 

Remark 4.2 Under a weaker assumption on the regularity of the solution, i.e. there 
exists an integer TOO > 4 such that (j) G Xj/ G C{[0,T];Hp°^^ ('^)) 

ncl[0,r];//7^‘(f2))nc2([0,r];//7^^(f2)) and 

ll^llL^florl'/r^o+S) + ^ II'^'^IIl”([o rl'/r^o+S) ^ I 

2 (4.65) 

II'^IIl°°([0,7’];//"'0+1) + II *^2'V^llz.“([0,7’];//"'0-3j ^ 1 , 

we can establish an //'-norm estimate of the MTI-FP method in Theorem 4.1 by a 
very similar proof with all the //^-norm in the proof being changed into //'-norm. 
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In ID case, the //'-norm estimate of the MTI-FP method holds without any stability 
(or CFL) condition. However in 2D/3D, due to the use of the inverse inequality to 
provide the bound in /“-norm of the numerical solution [6-8], we need impose the 
technical condition 

T<Pdih), With = ^ = (4.66) 

Of course, if the solution of the KGS is smooth enough, we can always adapt the 
//^-norm estimate in Theorem 4.1 and thus there is no need to assume the stability 
(or CFL) condition (4.66). 


5 Numerical results 

In this section, we report numerical results to demonstrate the uniform convergence 
of the MTI-FP method for e G (0,1] and apply it to numerically study convergence 
rates of the KGS equations to its limiting models ( 1 .9) and (1 .1 1) in the nonrelativistic 
limit regime. In order to do so, we take d = I and /i = A = lin(1.3) and choose the 
initial data in (1.5) with c/ = 1 as 

1 -f / / X 1 ^ 1 2 

roW = ^-sech ( y 1 , (j)o{x) =, 0i (x) = , x G K. (5.1) 


5.1 Accuracy test 

The problem (1.3) with d =l and (5.1) is solved on a bounded interval Q = [—32,32], 
i.e. b = —a = 32, which is large enough to guarantee that the periodic boundary con¬ 
dition does not introduce a significant aliasing error relative to the original problem. 
To quantify the error, we introduce two error functions: 

e^;''(f = tn) := := max^ , 

= tn) ■■= Wwi-dn) - Wi\\h2 , e5;'L(f) := ^max^ {^w!e{t)} . 

Since the analytical solution to this problem is not available, so the ‘reference’ solu¬ 
tion here is obtained numerically by the MTI-FP method (3.15)-(3.17) with very fine 
mesh h = 1/32 and time step T = 5 X 10 Tab. 5.1 and Tab. 5.2 show the spatial 
error of the MTI-FP method for \j/ and 0, respectively, at f = 1 under different e and 
h with a very small time step T = 5 x 10^® such that the discretization error in time 
is negligible. Tab. 5.3 and Tab. 5.4 show the temporal error of the MTI-FP method 
for \j/ and 0, respectively, at f = 1 under different £ and T with a small mesh size 
h= 1/16 such that the discretization error in space is negligible. 

From Tabs. 5.1-5.4 and extensive additional results not shown here for brevity, 
we can draw the following observations: 

(i) The MTI-FP method is spectrally accurate in space, which is uniformly for 
0 < £ < 1 (cf. Tabs. 5.1&5.2). The approximation in 0 is more accurate than y/. 
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Table 5.1: 

Spatial error analysis 

II 

1) with T = 5 X 10 

for different £ and h. 

II 

14-S- 

QJ 

ho = 1 

fto/2 

/io/4 

(io/8 

/!o/16 

eo = 1 

6.90E-1 

1.39E-1 

1.70E-3 

2.44E-7 

8.96E-11 

eo/2 

7.57E-! 

1.40E-1 

1.70E-3 

2.40E-7 

8.37E-11 

eo/22 

6.37E-1 

1.45E-1 

1.70E-3 

2.40E-7 

9.19E-11 

eo/23 

6.46E-1 

1.30E-1 

1.90E-3 

2.40E-7 

8.89E-11 

£t)/2'‘ 

6.47E-1 

1.31E-1 

1.80E-3 

2.38E-7 

9.26E-11 

eo/25 

6.47E-1 

1.31E-1 

1.80E-3 

2.38E-7 

9.66E-11 

eo/2'’ 

6.46E-1 

1.27E-1 

1.70E-3 

2.52E-7 

8.76E-11 

eo/2’ 

6.46E-1 

1.31E-1 

1.80E-3 

2.38E-7 

9.08E-11 

eo /28 

6.46E-1 

1.27E-1 

1.70E-3 

2.52E-7 

9.08E-11 

eo/2‘’ 

6.46E-1 

1.31E-1 

1.80E-3 

2.38E-7 

1.02E-10 

eo/2“ 

6.46E-1 

1.31E-1 

1.80E-3 

2.38E-7 

9.07E-11 

eo/2'3 

6.46E-1 

1.31E-1 

1.80E-3 

2.38E-7 

8.61E-11 


Table 5.2: Spatial error 

analysis: ej’* (f = 1) 

with T = 

5x10 ® for different £ and h. 


ho = 1 

ho/2 

/!o/4 

ho/S 

/!o/16 

eo = 1 

4.74E-2 

1.60E-3 

3.64E-6 

3.57E-10 

9.23E-11 

eo/2 

8.29E-2 

1.90E-3 

1.23E-5 

3.92E-10 

6.18E-11 

eo/2’ 

2.32E-1 

1.12E-2 

4.61E-5 

1.46E-9 

5.24E-11 

eo/23 

3.49E-1 

1.30E-3 

9.31E-5 

5.33E-10 

4.58E-11 

et)/2^ 

2.64E-1 

1.20E-3 

3.35E-5 

5.16E-9 

4.49E-11 

eo/2^ 

2.71E-1 

2.09E-4 

6.92E-6 

8.76E-10 

5.49E-1! 

eo/2'’ 

1.13E-1 

5.27E-4 

2.56E-6 

1.73E-10 

6.04E-1! 

eo/2’ 

3.99E-1 

4.29E-4 

6.28E-7 

4.23E-11 

6.20E-11 

eo/2* 

6.74E-2 

5.90E-4 

7.20E-8 

2.38E-11 

6.20E-11 

eo/2‘’ 

3.37E-1 

5.13E-5 

2.78E-8 

1.86E-11 

6.51E-11 

eo/2“ 

3.09E-1 

7.75E-4 

2.98E-9 

1.73E-11 

6.19E-1! 

eo/2‘’ 

1.87E-1 

3.64E-4 

3.33E-9 

1.72E-11 

6.13E-1! 


(ii) The MTI-FP method converges uniformly in time for e € (0,1] with lin¬ 
ear/quadratic al convergence rate for y//^, respectively (cf. last row in Tabs. 5.3(&5.4). 
In addition, for each fixed £ = eo > 0, when T is small enough, it converges quadrat- 
ically in time (cf. each row in the upper triangle of Tabs. 5.3&5.4); for fixed T, when 
0 < £ ^ 1, it converges quadratically in terms of £ (cf. each column in the lower 
triangle of Tabs. 5.3&5.4). 

(iii) The practical temporal uniform convergence order of the MTI-FP in (j) is 
better than the theoretical error estimate result. In fact, the temporal local truncation 
error of the MTI-FP in 0 is at 0{z^), but due to the nonlinear coupling with iff, the 
theoretical error bound drops down to uniform linear convergence rate. The rigorous 
mathematical analysis is still on-going. 

(iv) In the nonrelativistic limit regime, i.e. 0 < £ <C I, the mesh strategy (or £- 
scalability) of the MTI-FP method is T = (9(1) and h = (9(1). 
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Table 5.3: Temporal 
different e and T. 

error analysis: 

= 1) and 

=1) with h = 

1/16 for 


To = 0.2 

To/22 

To/2^ 

To/2« 

To/2** 

to/2‘0 

to/2>2 

£0 = 1 
rate 

1.31E-1 

1.60E-2 

1.52 

l.OOE-3 

2.00 

6.20E-5 

2.01 

3.84E-6 

2.01 

2.35E-7 

2.02 

1.77E-8 

1.87 

£o/2 

rate 

1.59E-1 

2.05E-2 

1.48 

1.30E-3 

1.99 

7.70E-5 

2.04 

4.76E-6 

2.01 

2.92E-7 

2.01 

1.37E-8 

2.20 

rate 

1.20E-1 

1.48E-2 

1.35 

1.20E-3 

1.81 

7.65E-5 

1.99 

4.79E-6 

2.00 

2.95E-7 

2.01 

1.38E-8 

2.20 

W2^ 

rate 

2.23E-2 

5.30E-3 

1.04 

5.19E-4 

1.68 

3.20E-5 

2.00 

2.00E-6 

2.00 

1.23E-7 

2.01 

5.77E-9 

2.20 

eo/2‘^ 

rate 

9.50E-3 

3.50E-3 

0.72 

2.20E-3 

0.33 

1.32E-4 

2.03 

8.06E-6 

2.02 

4.94E-7 

2.01 

2.31E-8 

2.20 

eo/2^ 

rate 

1.40E-3 

l.OOE-3 

0.24 

6.34E-4 

0.33 

4.67E-4 

0.22 

3.35E-5 

1.91 

2.02E-6 

2.02 

9.41E-8 

2.20 

^2^ 

rate 

3.82E-4 

5.65E-5 

1.38 

1.43E-4 

-0.67 

1.51E-4 

-0.18 

1.16E-4 

0.19 

8.29E-6 

1.90 

3.80E-7 

2.22 

^2^ 

rate 

5.64E-5 

2.30E-5 

0.65 

1.28E-5 

0.42 

3.61E-5 

-0.75 

3.83E-5 

-0.04 

2.89E-5 

0.20 

1.58E-6 

2.10 

W2® 

rate 

2.10E-5 

5.60E-6 

0.95 

6.41E-6 

-0.10 

2.61E-6 

0.65 

l.lOE-5 

-1.00 

1.12E-5 

-0.01 

5.70E-6 

0.49 

rate 

5.82E-6 

4.05E-6 

0.26 

2.40E-6 

0.38 

2.00E-6 

0.13 

2.41E-6 

-0.13 

8.49E-7 

0.75 

2.50E-7 

0.88 

£ 0 / 2 " 

rate 

2.09E-7 

2.19E-7 

-0.04 

1.67E-7 

0.20 

1.58E-7 

0.04 

8.48E-8 

0.45 

1.43E-7 

-0.38 

4.63E-8 

0.81 

rate 

1.78E-8 

8.93E-9 

0.50 

5.82E-9 

0.31 

5.68E-9 

0.02 

7.48E-9 

-0.20 

1.92E-9 

0.98 

9.46E-9 

-1.15 

e5;[L(r = 1) 

rate 

1.59E-1 

2.05E-2 

1.48 

2.20E-3 

1.61 

4.67E-4 

1.11 

1.16E-4 

1.00 

2.89E-5 

1.00 

5.70E-6 

1.17 


5.2 Convergence rates of KGS to its limiting models 


Let \j/ and (j) be the solution of the problem (1.3) with d = ji = X = I and (5.1), 
which is obtained numerically on a bounded interval Q = [—512,512] by the MTl-FP 
method (3.15)-(3.17) with very fine mesh h = l/16 and time step T = 10^^. Similarly, 
let t/sw and Zsw be the solution of the problem (1.9) with d = jX = X = 1 and (1.10) 
and (5.1), which is obtained numerically on the bounded interval Q by the EWI-FP 
method [7] with very fine mesh h = 1/16 and time step T = 10^^; and let Xffs and Zs 
be the solution of the problem (1.11) with d = jJ. = X = 1 and (1.12) and (5.1), which 
is obtained numerically on the bounded interval Q by the TSFP method [5] with very 
fine mesh li= 1/16 and time step T = 10^^. Denote 


0sw(jc,f) 




and define the error functions as 


t7sw(0 •—ll0(’iO 0sw(‘,OII//'+ IIV4w(‘,0ll/ri 1 
77s(f) := ||0(-,f) - (/)s(-,f)ll//i + II - rs(-,f)ll//i ■ 


(5.2) 
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Table 5.4: 

Temporal error analysis: el’g(f 

= 1) and 


with h = 

1/16 for 

different £ 

and T. 








To = 0.2 

To/22 

ro/2^ 

To/2'5 

To/2*' 

to/2‘0 

To/2‘2 

£ 0 = 1 

3.74E-2 

1.50E-3 

3.82E-5 

1.78E-6 

1.09E-7 

6.73E-9 

3.32E-10 

rate 

— 

2.32 

2.65 

2.21 

2.01 

2.01 

2.17 

£ 0/2 

9.32E-2 

4.80E-3 

2.29E-4 

1.35E-5 

8.32E-7 

5.11E-8 

2.39E-9 

rate 

— 

2.14 

2.19 

2.04 

2.01 

2.01 

2.20 

V2^ 

9.65E-2 

1.26E-2 

7.09E-4 

4.23E-5 

2.61E-6 

1.60E-7 

7.5 IE-9 

rate 

— 

1.47 

2.07 

2.03 

2.01 

2.01 

2.20 

V2^ 

2.33E-2 

4.10E-3 

3.02E-4 

1.95E-5 

1.21E-6 

7.47E-8 

3.50E-9 

rate 

— 

1.25 

1.88 

1.98 

2.00 

2.01 

2.20 

eo/2* 

8.70E-3 

3.30E-3 

9.09E-5 

6.00E-6 

3.82E-7 

2.36E-8 

l.llE-9 

rate 

— 

0.70 

2.59 

1.96 

1.99 

2.00 

2.20 

eo/2^ 

9.59E-4 

8.29E-4 

2.42E-4 

3.73E-6 

2.89E-7 

1.81E-8 

8.53E-10 

rate 

— 

0.11 

0.89 

3.01 

1.85 

2.00 

2.20 


4.71E-4 

1.19E-4 

1.33E-5 

3.80E-6 

2.13E-7 

1.67E-8 

8.01E-10 

rate 

— 

0.99 

1.58 

0.90 

2.07 

1.84 

2.19 


8.27E-5 

2.88E-5 

2.49E-6 

1.76E-7 

4.95E-8 

1.30E-8 

7.86E-10 

rate 

— 

0.76 

1.76 

1.91 

0.92 

0.96 

2.02 

W2® 

3.33E-5 

5.69E-6 

3.28E-7 

2.18E-8 

3.37E-9 

1.70E-9 

6.57E-10 

rate 

— 

1.28 

2.06 

1.95 

1.35 

0.49 

0.69 

Wz’ 

8.92E-6 

3.98E-6 

1.04E-7 

5.46E-9 

3.66E-10 

5.75E-11 

2.04E-11 

rate 

— 

0.58 

2.63 

2.12 

1.95 

1.34 

0.75 

£ 0 / 2 " 

3.05E-7 

1.59E-7 

5.64E-9 

1.32E-9 

7.21E-11 

7.09E-11 

7.43E-11 

rate 

— 

0.47 

2.41 

1.05 

2.10 

0.01 

-0.03 


2.30E-8 

7.26E-9 

2.21E-9 

7.77E-10 

7.74E-10 

7.74E-10 

7.74E-10 

rate 

— 

0.83 

0.86 

0.76 

0.01 

0.00 

0.00 


9.65E-2 

1.26E-2 

7.09E-4 

4.23E-5 

2.61E-6 

1.60E-7 

7.51E-9 

rate 

— 

1.46 

2.00 

2.00 

2.00 

2.00 

2.00 

Fig. 5. 

1 depicts T], 

iw(0 and 

Tjsit) for different £ 

. From Fig. 

5.1 and additional 


similar results not shown here for brevity, when £ 0, we can observe numerically 

the solution of the KGS equations (1.3) converges quadratically and uniformly in 
time to that of its limiting model (1.9) as 

Wvi'if) ~ ’/sw{‘,0ll//i(i2) + II0(‘!O ~ 0sw(‘,OII/ri(t2) — t >0, (5.3) 

and respectively, to its limiting model ( 1.11 ) as 

IIV^(-,O-rs(-,f)ll//>(r2) + ll0(-,O-0s(-,OII//i(i2) < {Ci+C2T)e\ 0<t<T, 

(5.4) 

where Co, Ci and C 2 > 0 are positive constants which are independent of £ and T. 
Rigorous mathematical justification of these numerical observations is on-going. 


6 Conclusion 

A multiscale time integrator Fourier pseudospectral (MTI-FP) method was proposed 
and analyzed for solving the Klein-Gordon-Schrodinger (KGS) equations with a di- 
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t 




Fig. 5.1: Time evolution of TjswlO and 77s(0 for different e. 


mensionless parameter 0 < e < 1 which is inversely proportional to the speed of 
light. The key ideas for designing the MTI-FP method are based on (i) carrying out 
a multiscale decomposition by frequency at each time step with proper choice of 
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transmission conditions between adjacent time intervals, and (ii) adapting the Fourier 
spectral method for spatial discretization and the EWI for integrating second-order 
highly oscillating ODEs. Rigorous error bounds for the MTI-EP method were estab¬ 
lished, which imply that the MTl-EP method converges uniformly and optimally in 
space with spectral convergence rate, and uniformly in time with linear convergence 
rate for £ G (0,1] and optimally with quadratic convergence rate in the regimes when 
either £ = (9(l)orO<£<T. Numerical results confirmed these error bounds and 
showed that the MTl-EP method offers compelling advantages over classical methods 
in the nonrelativistic regime 0 < £ <C 1. Numerical results suggest that the solution of 
the KGS equations converges quadratically to that of its limiting models when £ —^ 0. 
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Appendix A. Explicit formulas for the coefficients used in (3.14) 

Define the functions 



(A.l) 


Taking s = T in (A.l), after a detailed computation, we have 
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